# chooseBioCmirror()
# BiocManager::install('GreenleafLab/ArchR', ref = 'master')
library(ArchR)
# installExtraPackages()
set.seed(1)

inputFiles <- getTutorialData("Hematopoiesis")
inputFiles

# Ref genome data is needed
# need to be in accord with fragment file
addArchRGenome("hg19")

# create arrow files ---------
# take ~5 min
ArrowFiles <- inputFiles |>
  createArrowFiles(
  minTSS = 4, #Dont set this too high because you can always increase later
  minFrags = 1000, 
  addTileMat = TRUE,
  addGeneScoreMat = TRUE
)

ArrowFiles

# inferring doublets ---------
doubScores <- ArrowFiles |>
  addDoubletScores(
  k = 10, #Refers to how many cells near a "pseudo-doublet" to count.
  knnMethod = "UMAP", #Refers to the embedding to use for nearest neighbor search.
  LSIMethod = 1
)

# ArchRProject --------
proj <- ArrowFiles |>
  ArchRProject(
  outputDirectory = "HemeTutorial",
  copyArrows = TRUE #This is recommened so that you maintain an unaltered copy for later usage.
)

getAvailableMatrices(proj)

proj <- filterDoublets(proj)

# dimreduc & clustering ------
proj <- addIterativeLSI(proj) |>
  addClusters() |>
  addUMAP()

p1 <- proj |>
  plotEmbedding(colorBy = "cellColData", name = "Sample", embedding = "UMAP")

p2 <- proj |>
  plotEmbedding(colorBy = "cellColData", name = "Clusters", embedding = "UMAP")

ggAlignPlots(p1, p2, type = "h")

plotPDF(p1,p2, name = "Plot-UMAP-Sample-Clusters.pdf",
        ArchRProj = proj, addDOC = FALSE, width = 5, height = 5)

## native plot func failed
## try plot manually ------
tidy.arc <- proj@cellColData |>
  as_tibble(rownames = '.cell')

tidy.arc <- getEmbedding(proj) |>
  as_tibble(rownames = '.cell') |>
  left_join(tidy.arc)

tidy.arc |>
  ggplot(aes(-`IterativeLSI#UMAP_Dimension_1`, -`IterativeLSI#UMAP_Dimension_2`,
             color = Sample)) +
  geom_point()

tidy.arc |>
  ggplot(aes(`IterativeLSI#UMAP_Dimension_1`, `IterativeLSI#UMAP_Dimension_2`,
             color = Clusters)) +
  geom_point() +
  scale_color_brewer(palette = 'Paired')

# assign clusters with genetic scores ------
proj <- addImputeWeights(proj)

markerGenes  <- c(
  "CD34",  #Early Progenitor
  "GATA1", #Erythroid
  "PAX5", "MS4A1", "MME", #B-Cell Trajectory
  "CD14", "MPO", #Monocytes
  "CD3D", "CD8A"#TCells
)

p <- plotEmbedding(
  ArchRProj = proj, 
  colorBy = "GeneScoreMatrix", 
  name = markerGenes, 
  embedding = "UMAP",
  imputeWeights = getImputeWeights(proj)
)

#Rearrange for grid plotting
p2 <- lapply(p, function(x){
  x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
      axis.text.x=element_blank(), 
      axis.ticks.x=element_blank(), 
      axis.text.y=element_blank(), 
      axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

# plot genome browser tracks ----------
p <- plotBrowserTrack(
  ArchRProj = proj, 
  groupBy = "Clusters", 
  geneSymbol = markerGenes, 
  upstream = 50000,
  downstream = 50000
)

grid::grid.newpage()
grid::grid.draw(p$CD14)

## save multi-page pdf ------
plotPDF(plotList = p, 
        name = "Plot-Tracks-Marker-Genes.pdf", 
        ArchRProj = proj, 
        addDOC = FALSE, width = 5, height = 5)

# shiny genome browser ------
ArchRBrowser(ArchRProj = proj)

proj <- saveArchRProject(ArchRProj = proj)
proj <- loadArchRProject(path = "HemeTutorial")
